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Valery WebeiQ 

Department of Chemistry, University of Fribourg, 1700 Fribourg, Switzerland. 



o 
o 

(N 
o 

Q 

(N 



C/3 



Anders M. N. Niklasson and Matt Challacombe 
Los Alamos National Laboratory, Theoretical Division, Los Alamos 87545, New Mexico, 

(Dated: February 2, 2008) 



USA. 



X3 

o 



> 

00 

o 



o 



X 



Perturbed projection for linear scaling solution of the coupled-perturbed self-consistent-field equa- 
tions [Weber, Niklasson and Challacombe, Phys. Rev. Lett. 92, 193002 (2004)] is extended to the 
computation of higher order static response properties. Although generally applicable, perturbed 
projection is developed here in the context of the self-consistent first and second electric hyperpolar- 
izabilities of three dimensional water clusters at the Hartree-Fock level of theory. Non-orthogonal, 
density matrix analogues of Wigner's 2n + 1 rule are given up to fourth order. Linear scaling 
and locality of the higher order response densities under perturbation by a global electric field are 
demonstrated. 

PACS numbers; 02.70.-c, 71.15.Dx,31.15.Ar,31.15.Md,31.15.Ne,33.15.Kr,36.40.Cg 



INTRODUCTION 



First principles electronic structure theory has tradi- 
tionally been limited to the study of small systems with 
a limited number of nonequivalent atoms. Despite the 
tremendous increase in computational power of digital 
computers this has remained the case, until the advent 
of reduced complexity algorithms over the last decade 

1^1 IE Q Hi 111 ' In the best case, these reduced complex- 
ity algorithms scale only linearly with system size, TV, al- 
lowing simulation capabilities to keep pace with hardware 
improvements. These linear scaling algorithms exploit 
the quantum locality (or nearsightedness) of non-metallic 
systems, manifested in the approximate exponential de- 
cay of density matrix elements with atom-atom separa- 
tion through the effective use of sparse matrix methods. 
For small systems, linear scaling methods may be ineffi- 
cient due to overhead. However, for large, complex sys- 
tems these methods hold the promise of major impact 
across materials science, chemistry and biology. 

So far, a majority of work in linear scaling electronic 
structure theory has focused on methods and calculations 
involving the ground state, with little attention devoted 
to the problem of response properties. The calculation 
of static response within Hartree-Fock or Density Func- 
tional Theory may be obtained through solution of the 
Coupled-Perturbed Self-Consistent-Field (CPSCF) equa- 
tions, which yield properties such as the electric polariz- 
ability and hyperpolarizability 0, j the Born-effective 
charge, the nuclear magnetic shielding tensor 0, indirect 
spin-spin coupling constant |lOl llll |. geometric deriva- 
tives (i.e. higher order analytic force constants) and 
polarizability derivatives such as the Raman intensity 
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[l^l Il3|. to name but a few. 

Conventional ap proaches to solution of the CPSCF 
equations |2lla,ll3 based on perturbation of the wave 
function, requiring an 7V"^-scaling eigensolve which may 
need to be followed by an 0{N^) transformation of two- 
electron integrals, depending on the method. In addition 
to the formal scaling of these conventional methods, they 
do not admit exploitation of quantum locality through 
the effective use of sparse matrix algebra. More recently, 
schemes with the potential for reduced complexity have 
been been put forward. Ochsenfeld and Head-Gordon 
proposed a scheme based on the Li-Nunes-Vanderbilt 
density matrix minimization [Tel. Later, Larsen et al. 
[itI i proposed iterative solution of the CPSCF equations 
involving equations derived from unitary operations and 
approximations to the matrix exponential. In both of 
these approaches, a linear system of equations containing 
commutation relations is obtained, which implicitly de- 
termines the response function. However, the method of 
solution for these equations is not discussed, and compu- 
tational results are not presented. Recently though, with 
an apparent reformulation of Ref. JjGj, Ochsenfeld, Kuss- 
mann and Koziol |l8l | have achieved near linear scaling 
computation of NMR chemical shifts for one-dimensional 
alcanes at the GIAO-HF/6-31G* level of theory, but 
likewise provide no details on their method of solution. 
These implicit commutation relations are Sylvester-like 
and may be solved with a number of approaches |1^] , the 
particulars of which arc of interest. 

In contrast. Perturbed Projection is a recently 
developed alternative for Ai"-scaling solution of the CP- 
SCF equations that is simple and explicit. Based on 
a recently developed density matrix perturbation the- 
ory [ij. Perturbed Projection exploits the explicit re- 
lationship between the density matrix V and the effec- 
tive Hamiltonian or Fockian J- via spectral projection; 
V = 0{jll — T), wherein is the Heaviside step func- 
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tion (spectral projector) and the chemical potential /2 
determines occupied states via Aufbau filling. Spec- 
tral pro ject ion can be carried out in a number of ways 
miillillilSliSlillii- Of special interest here 
are recursive polynomial expansions of the projector, in- 
cluding the second order trace correcting (TC2) and 
fourth order trace resetting (TRS4) purification algo- 
rithms. These new methods (TC2 and TRS4) have con- 
vergence properties that depend only weakly on the band 
gap, do not require knowledge of the chemical potential 
and perform well for all occupation to state ratios. Per- 
haps most important to the current contribution, these 
methods converge rapidly to smooth, monotone projec- 
tors. 

Prior to Ref. Perturbed Projection demonstrated 
linear scaling in computation of the first electric polar- 
izability for three-dimensional water clusters with the 
Hartree-Fock model [l^. Also, in a preceding paper, 
we outlined a non-orthogonal density matrix perturba- 
tion theory [sOl for response to a change in basis (i.e. as 
occurs in the evaluation of higher order geometric energy 
derivatives In this article, the Perturbed Projec- 

tion method is extended to higher orders in the electric 
polarizability, up to fourth order in the total energy. 

This paper is organized as follows: First we describe 
the perturbation expansion and the computation of re- 
sponse properties through solution of the CPSCF equa- 
tions. Then we present extension of Perturbed Projection 
through higher orders and the computation of properties 
using a density matrix analogue of Wigner's 2n + 1 rule. 
Next, we present several examples of calculated higher 
order response properties. We show the saturation of 
hyperpolarizabilities up to fourth order (i.e. up to the 
second hyperpolarizability 7) for a series water chains. 
We also demonstrate linear scaling complexity for the 
solution of the higher order CPSCF equations and an 
approximate exponential decay in elements of higher or- 
der response functions for 3D water clusters. Finally, we 
discuss these results and present our conclusions. 

THE COUPLED PERTURBED 
SELF-CONSISTENT-FIELD EQUATIONS 

The Coupled-Perturbed Self-Consistent-Field (CP- 
SCF) equations yield static response functions and prop- 
erties in models including both the Hartrcc-Fock (HF) 
and Density Functional Theory (DFT). In the follow- 
ing we develop Perturbed Projection for solution of the 
CPSCF equations in the framework of polarization and 
Hartree-Fock theory. In many cases, the extension of 
Perturbed Projection to the computation of other static 
perturbations is straightforward. In the case of model 
chemistries that involve DFT, an extra programming ef- 
fort is required [^H^- Also, in the case of properties 
such as the NMR chemical shift and geometric derivatives 



(force constants) , perturbation of the non-orthogonal ba- 
sis requires additional considerations that we have de- 
tailed in a preceding paper [sot- 

Notation 

Superscripts and subscripts refer to perturbation order 
and self-consistent cycle count respectively. The symbols 
T>,J-,... are matrices in an orthogonal representation, 
while D, F, . . . are the corresponding matrices in a non- 
orthogonal basis. The transformation between orthogo- 
nal and non-orthogonal representations is carried out in 
0{N) using congruence transformations provided 
by the approximate inverse (AINV) algorithm for com- 
puting sparse approximate inverse Cholesky factors with 
a computational complexity scaling linearly with the sys- 
tem size Umillll- 

Response expansions 

Within HF theory, the total electronic energy Etot of 
a molecule in a static electric field £ is 

Etot{£) = Tr{D ■ (/jO + ^l£)} + ]^Tr{D ■ {J[D] + K[D])} 

= Tt{D ■ F[D]} - \Tr{D ■ {J[D] + K[D])}, 

(1) 

where D = D[£] is the density matrix in the electric field 
£, is the core Hamiltonian, /i is the dipole moment 
matrix, J[D] is the Coulomb matrix, K[D\ the exact HF 
exchange matrix and 

F = F[£] = h'^ + li£ + J[D{£)]+K[D{£)] (2) 

is the Fockian. The total energy of a molecule in a ho- 
mogeneous electric field may be developed in a Taylor 
expansion series around £ — Q as, 

£;tot(f)=Stot(0)-^Maf" 

a 

-^Y.^a,£''£' 

ah 

-^^/Safccf^f'f^ (3) 
ahc 

abed 

+ ..., 

where Uab is the polarizability, and (3abc and ^abcd are the 
first and second hyperpolarizabilities, respectively, ^la is 
the dipole moment, and f ° is the electric field in direction 
a. The polarizability aab is the second order response of 
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the total energy with respect to variation in the electric 
field while the higher derivatives, Pabc and 7ahi-rfj_Kive rise 
to the first and second hyperpolarizabilities jll 1^ where 



Ctab 



Pabc = 



= -2Tr[Dyb], 



(4a) 



£^0 



^abcd 



d^Etot 



= -4Tr[D''''fic], (4b) 



£=0 



dE-^d£''dE''dE'^ 



= -l2Tr[D''''^id]. (4c) 



£=0 



Here - denotes a density matrix derivative with re- 
spect to a field in directions a . . . at £ = 0. The density 
matrix derivative or "response function" is given by 



pa... ^ 



—e{pi-T{E)) 



d£ 



(5) 



£=0 



The Fockian may also be expanded order by order in the 
pertm'bation to yield 



T{£) 



a 

2! ^ 



ah 



(6) 



ihc 



where stands for dT{E)/dE'', T'''' = d^T{E) / dE'^dE'' , 
and so on for the higher order terms. A similar expansion 
also holds for the density matrix !?(£). 

Within HF theory, the unperturbed Fockian in the 
non-orthogonal basis is 



F° + J{D°)+K{D°), 
while the first variation of the Fockian is 

F^ = ^ia + JiD") + KiD") 
and the higher terms are given by 

pab... ^ j[D''^ -)+K{D'''' -). 



(7) 



(8) 



(9) 



In computation of the unperturbed Fockian, the Coulomb 
matrix J may be computed in 0{N\^N) with the Quan- 
tum Chemical Tree Code (QCTC) ,3^] and the exchange 
matrix K computed in 0{N) with the C'(A^)-exchange 
(ONX) algorithm that exploits quantum locality of the 
density matrix D° j^^- Likewise the Fockian derivatives, 
F" ", may be computed with the same algorithms in lin- 
ear scaling time if elements of _D°- - manifest an approxi- 
mate exponential decay with atom-atom separation, sim- 
ilar to the decay properties of D^ . 

While the expansions above are given explicitly for 
Hartree-Fock Theory, similar expressions hold also for 
Kohn-Sham and hybrid HF/DFT, which involve varia- 
tion of the exchange-correlation matrix V^^' {D^, D'^, . . .) 

MM. 



Conditions for self-consistency 

The derivative density matrices and derivative Fock- 
ians depend on each other implicitly, and must be solved 
for self-consistently via the CPSCF equations. The nec- 
essary and sufficient criteria for convergence of the CP- 
SCF equations involve generalized self-consistence condi- 
tions 



= 0, 

[^=■",1?°] + [T^^V] =0, 



(10) 
(11) 



-f =0, 
3 3 

+ ^[:f^\v^] + ^[:f%v'^] ^^^^ 

3 3 
in addition to generalized idempotency-like constraints 



^ {X>",X'°}, 



'lb 

1 

3 



(14) 
(15) 

(16) 
(17) 



where the anti-commutator notation {A, B} — AB + BA 
has been used. 



SOLVING THE HIGHER ORDER CPSCF 
EQUATIONS WITH PERTURBED PROJECTION 

In solution of the CPSCF equations, it is first neces- 
sary to determine the ground state density matrix I?". 
This may be accomplished in 0(N) using a purification 
algorithm such as Niklasson's [221 second order trace cor- 
recting scheme (TC2) in conjunction with sparse atom- 
blocked linear algebra [23, |4l| . Linear scaling is achieved 
for insulating systems through the dropping (filtering) of 
atom-atom blocks with Frobenious norm below a numeri- 
cal threshold (r ^ 10^'* — 10^^). At SCF convergence the 
TC2 algorithm generates a polynomial sequence defin- 
ing the ground state projector, from which the derivative 
density matrices are directly obtained. 

Having solved the groundstate SCF equations, solution 
of the CPSCF equations commences with a guess at the 
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derivative densities, followed by computation of deriva- 
tive Fockians. At the r*-^ CPSCF cycle, the n"" order 
derivative Fockians are 



fia + J{D^) + K{D^), n=l 
J{D^-) + K{D^-), n>l. 



(18) 



After construction of the derivative Fockians, response 
functions through D!^_^i are computed, constituting one 
cycle in solution of the CPSCF. As described in Section 
B these response functions are obtained directly through 
variation of the occupied subspace projector, 



jja... _ 



-e{fLl - Tr{£)) 



(19) 



£=0 



which is accomplished via the Niklasson and Challa- 
combe density matrix perturbation theory [2^ . 

After a few CPSCF cycles, the approach to self- 
consistency may be accelerated with Weber and Daul's 
DDES algorithm 



E 

k—r~i 



(20) 



in which the Ck coefficients are chosen to minimize the 
n'^ order commutation relations, as in Eqs. (|ll() - (|13|l . 
The application of the DDIIS algorithm to acceleration 
of higher order CPSCF equations is developed further in 
Section B 

At self-consistency, the conditions given in Section b 
have been met, and it is then appropriate to compute re- 
sponse properties. In general, we can use the expectation 
value 



£;(')') = 2(7 - iy.Tr{D'^^-^^h^^'^), 



(21) 



where in the case of polarization we have the 
(hyper)polarizabilities Uat = — 2Tr[Z3°/ih], Pabc — 
-4Tr[D°Vc] or jabcd = -12Tr[i?'^''=/irf]. Alternatively, 
it is possible to construct density matrix analogues of 
Wigner's 2n + 1 rule, which allows the evaluation of re- 
sponse properties up to order 2n -I- 1 from response func- 
tions of order n. These analogues are given in Section b 
through 4"* order in the energy. 



Perturbed Projection 

Although a number of analytic, asymptotically dis- 
continuous representations exist for the Heaviside step 
function 9, direct representation (and variation) of these 
forms as in Eq. (|19() is problematic. Polynomial expan- 
sion of the step function is a alternative choice, but de- 
mands a very high order and can be costly. Specifically, 
polynomial expansion of 9 with a p'th order polynomial 
incurs a cost that is at best 0{y/p) 01- Polynomial 
expansion techniques, such as those based on the the 



Tchebychev polynomials, may also be plagued by Gibbs 
oscillations |4J|, which are high order ripples in the ap- 
proximate 9 due to incompleteness. Alternatively, recur- 
sive purification methods achieve high order representa- 
tion in 0{\ogp) j2^]. Also, purification methods (such 
as TC2 and TRS4) yield projectors that are smooth and 
strictly monotonic. 

Each Perturbed Projection sequence is based on a cor- 
responding purification scheme or generator, such as TC2 
22]. The Perturbed Projection sequence is obtained by 
collecting terms of the response order by order upon per- 
turbative expansion of its generator |21| . Perturbed Pro- 
jection provides explicit, recursive formulae for the con- 
struction of response functions, retaining the convergence 
properties, smoothness and montonicity of the generat- 
ing sequence. These explicit formulae stand in contrast 
to methods where the density matrix derivatives are im- 
plicitly defined as solutions to equations of Sylvester type 

'Ti,'!?,'!^. 

Sufficient to compute fourth order properties using the 
2n + 1 rule presented in Section b Perturbed Projection 
is outlined in the following for computation of the sec- 
ond order response function. The Perturbed Projection 
sequence is started with the A'q ", which are prepared 
from the Fockians JF", and T""^ by cornpressing their 
spectrum into the domain of convergence |23| using 



and 



A.r, 



J fy 



J~ rr 



J fr 



■pa. 



- T 



(22) 



(23) 



where Tmin and Tmax are upper and lower bounds to the 
eigenvalues of JF°. 

While Perturbed Projection can be formulated within 
any purification scheme, we focus here on the simple and 
efficient TC2 method 22]. Briefly, TC2 constructs a 
ground state projector through a series of trace correct- 
ing projections; when the trace is larger than iVg, is 
used to reduce the trace, and when the trace is less than 
Nc, 2x — is used to increase the trace. The resulting 
sequence of correcting projections yields a step at the 
correct chemical potential. Within this framework, the 
second order TC2 Perturbed Projection sequence is 



■yah 

■ya 

■yh 

or 

yah 

'^ i+i - 



{xtxf) 



2X-^-{{X-\Xf} 
2X^-{XlX^} 
2X--{X-,Xf} 
2Xf - [Xff 



Tr[<^0]>iVe (24) 



Tt[X^] < Ne. 



(25) 
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RESULTS 



As with the TC2 generator, the n'^ order response func- 
tions 



pa.. 



nl lim A""- 



(26) 



converge quadratically, reaching convergence when either 
the error e = |Tr[A'0] - A^e| + |Tr[A7]| + . . or the 
maximum element in the change SX"" ' = — 
falls below r, the atom-atom block drop tolerance de- 
scribed in Ref. [23|. As described more completely in 
Ref. 30], when the solution gets close to convergence, 
i.e. |Tr[A'0] - iVe| < e with e « lO^^ - lO^^^ we alternate 
the projection at each step, which protects the conver- 
gence under the incomplete sparse linear algebra. 



Derivative DIIS 

Direct inversion in the iterative subspace (DIIS), in- 
troduced some time ago by Pulay [iil l4fil|. accelerates 
convergence toward self-consistency. DIIS employs in- 
formation accumulated during preceding iterations to 
construct an effective Fockian !Fk at the k-ih SCF cy- 
cle, which minimizes the commutation error between the 
Fockian and the density matrix. The effective Fockian is 
then used instead of Tk to generate an improved density 
matrix. 

Recently, Weber and Daul have developed the Deriva- 
tive DIIS (DDIIS) scheme for accelerating convergence of 
the CPSCF equations Like DIIS, DDIIS is based on 
minimization of the Frobenious norm of an error matrix 



(27) 



where the e° ' 's are just the n-th order commutator re- 
lation of Eqs. H10I13|I (e.g. the first order error matrix 
is given by e° = [Tf ,V°] + [J^°,Pf]). The optimal co- 
efficients Ci are solutions to the quadratic programming 
problem 



inf 



- y 

2 ^ 

i.j—r — s 



(28) 



where elements of the B matrix are given by By = 
Tr[e^' (e"' )"^]. A working equation is then obtained 
through the associated Euler-Lagrange equation 



B 1 

1^ 



(29) 



where = (0, . . . , 0)-^ and 1 — (1, . . . , 1)"^ are vectors 
whose components are and 1 respectively and A is the 
Lagrange multiplier of the constraint ^Y^=n-ni^i ^ 
The set of linear equations is solved by inverting the left- 
hand side matrix. In the event of a singular or near 
singular matrix, the rank of Eq. H29(l is reduced by dis- 
carding the oldest entries (increasing s) until the linear 
system stabilizes. 



Density matrix formulation of Wigner's 2?! + 1 rule 

Wigner's 2n+l rule, traditionally predicated on deriva- 
tives of the wavefunction, yields order 2n -I- 1 in the en- 
ergy response from n-th order derivatives 0, 0| . A den- 
sity matrix analogue of Wigner's 2n + 1 rule for a sin- 
gle perturbation parameter was given to third order by 
McWeeny ji^ and up to fourth order by Niklasson and 
Challacombe in the orthogonal representation. For 
completeness, we present non-orthogonal generalizations 
up to fourth order, which may require less memory due 
to the less dense structure of non-orthogonal matrix in- 
termediates. 

The first and second order energy corrections are well 
known, corresponding simply to expectation values as in 
Eq. H21(l . Beyond second order, the 2n -f 1 rule offers 
a valuable alternative. The third order non-orthogonal 
contribution is 



^abc ^ 2 



P(a,&,c) 



Ti[[D'',D°]sSD''F'= 



(30) 



where P{a, b, c) stands for the permutation operator 
such that all permutations of a, b and c are made 
(e.g. P(a, b, c) generates the sum of all the six terms: 
(a,5, c), (a, c, 5), (6,a,c), (b,c,a), {c,a,b) and (c, 6, a)) 
and [A, B\s — ABB — BSA where S is the overlap ma- 
trix. Similarly, the fourth order non-orthogonal contri- 
bution is 



bed 



P(a,b,c,d) 



Tt[[D''\D°]sSD''F'^- 



(31) 



[D'',D°]sS{D'"'F'^ + D^'F"'^)] 



For the orthogonal case S ~ I and Z?", F'^, . . . are re- 
placed by JF'^, .... In most cases the complexity of 
these equations can be reduced by taking advantage of 
indicial symmetry; a, 5, c, and d represent the Cartesian 
directions x,y, z so that terms with indecies in the same 
direction simplify. For example, reduces to only 

one term requiring only 15 matrix multiplications. In 
the worst case, where all the directions are different, i.e. 
j^aabc j-Qj, ^j-^y other permutation of (a, a, 5, c)), the re- 
lation l|31|l reduces to include only 12 terms with 180 
matrix multiplications. Similar reductions of the com- 
putational cost also apply to Eq. H30|l . The number of 
matrix-matrix multiplies can be further reduced if one 
uses an orthogonal representation, but this typically in- 
volves matrix-matrix multiplies with more dense inter- 
mediates. 



RESULTS 

We have implemented these methods in the Mon- 
doSCF suite of linear scaling quantum chemistry pro- 
grams 'soj. The construction of the Fockian and deriva- 
tive Fockian was carried out using the linear scaling 
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RESULTS 



TABLE L The longitudinal polarizability, a^z, for water 
chains at the RHF/6-31G level of theory, computed with 
MONDOSCF using GOOD and TIGHT numerical thresholds, and 
also with the GAMESS quantum chemistry package |49|| . 



A^HaO 


GAMESS 


GOOD 


TIGHT 


1 


5.8136 


5.813620 


5.813588 


2 


6.3448 


6.345037 


6.344822 


3 


6.5844 


6.584658 


6.584435 


4 


6.7276 


6.727905 


6.727672 


5 


6.8226 


6.823290 


6.822857 


10 


7.0308 


7.031056 


7.030858 


15 


7.1047 


7.104904 


7.104770 


20 


7.1424 


7.142580 


7.142422 



TABLE IL The longitudinal first hyperpolarizability, Pzzz, 
for water chains at the RHF/6-31G level of theory, computed 
with MONDOSCF using GOOD and TIGHT numerical thresholds, 
and also with the GAMESS quantum chemistry package 0|. 



A^HaO 


GAMESS 


GOOD 


GOOD" 


TIGHT 


TIGHT" 


1 


-30.6125 


-30.611029 


-30.612627 


-30.612163 


-30.612256 


2 


-29.5444 


-29.547427 


-29.548504 


-29.544907 


-29.544994 


3 


-25.3696 


-25.372208 


-25.373615 


-25.370297 


-25.370381 


4 


-22.1411 


-22.143436 


-22.145040 


-22.141494 


-22.141603 


5 


-19.8925 


-19.902088 


-19.904449 


-19.896462 


-19.897141 


10 


-14.8063 


-14.807075 


-29.617990 


-14.806973 


-14.807119 


15 


-12.9713 


-12.969238 


-12.972227 


-12.971940 


-12.972124 


20 


-12.0334 


-12.028709 


-12.033633 


-12.034014 


-12.034238 



"The density matrix based 2n+l rule has been used. 



QCTC method for computation of the Coulomb matrix 
[i^ Isij and the ONX algorithm [s^ [s^l for computa- 
tion of the Hartree-Fock exchange matrix. The CPSCF 
equations were solved in an entirely orthogonal represen- 
tation, with higher order properties evaluated using the 
2n + 1 rule with the non-orthogonal formulae given by 
Eqs. H30I31|I . Two different levels of numerical accuracy 
have been used, GOOD and TIGHT. Thresholds that de- 
fine the GOOD accuracy level include a matrix threshold 
r = 10~^. as well as other numerical thresholds detailed 
in Ref. |5l|, which deliver 6 digits of relative accuracy 
in the total energy. The TIGHT option involves the ma- 
trix threshold t = 10~^ and delivers 8 digits of relative 
accuracy in the total energy. 

Calculations were carried out on a single Intel Xeon 
2.4GHz processor running RcdHat Linux 8.0 and exe- 
cutables compiled with Portland Group Fortran Com- 
piler pgf90 4.0-2 

Convergence of the CPSCF equations for the water 
systems described in the following are typically achieved 
in about 10 cycles, independent of cluster size, basis set, 
matrix threshold or order of the response calculated. 

All results are reported in atomic units. Also, un- 



TABLE IIL The longitudinal second hyperpolarizability, 
^zzzz, for water chains at the RHF/6-31G level of theory, 
computed with MondoSCF using GOOD and TIGHT numerical 
thresholds^ and also with the GAMESS quantum chemistry 
package 49] . 

Nh^O GAMESS GOOD GOOD" TIGHT TIGHT" 

1 330.5753 330.54375 330.54319 330.57193 330.5724 

2 820.1398 820.19231 820.19667 820.14775 820.1493 

3 1008.5656 1008.5855 1008.6073 1008.5752 1008.5765 

4 1103.4813 1103.5053 1103.5280 1103.4883 1103.4902 

5 1168.9563 1169.2104 1169.2531 1169.0630 1169.0754 
10 1324.2906 1325.1208 1324.6321 1324.2975 1324.2999 
15 1381.8657 1383.0802 1382.2322 1381.8758 1381.8738 
20 1411.4264 1414.4092 1411.8528 1411.4410 1411.4292 

"The density matrix based 2n-|-l rule has been used. 



FIG. 1: Total CPU time of the fifth CPSCF iteration of fourth 
order for the water cluster sequence with the 6-31G and 6- 
31G** basis sets and the GOOD and TIGHT numerical thresholds 
(see text) controlfing numerical precision of the result. The 
lines are fits to the last three and four points, respectively. 




number of water molecules 



less otherwise noted, all timings and values have been 
obtained by computing the n"^ order response function 
and evaluation with the n + I rule (expectation value), 
Eq. 



One dimensional water chains 

Perturbed Projection has been used to compute the 
(hyper)polarizabilities azz, Pzzz and ^zzzz of linear wa- 
ter chains up to (H20)2o- These calculations have been 
carried out with MoNDoSCF at the RHF/6-31G level of 
theory using both the GOOD and TIGHT thresholding pa- 
rameters, as well as with the conventional algorithms im- 
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RESULTS 



FIG. 2; QCTC CPU time of the fifth CPSCF iteration of 
fourth order for the water cluster sequence with the 6-31G 
and 6-31G** basis sets and the GOOD and TIGHT numerical 
thresholds (see text) controlling numerical precision of the 
result. The lines are fits to the last three and four points, 
respectively. 



FIG. 4: ONX CPU time of the fifth CPSCF iteration of fourth 
order for the water cluster sequence with the 6-31G and 6- 
31G** basis sets and the GOOD and TIGHT numerical thresholds 
(see text) controlling numerical precision of the result. The 
lines are fits to the last three and four points, respectively. 



600 




number of water molecules 



FIG. 3: TC2 CPU time of the fifth CPSCF fteration of fourth 
order for the water cluster sequence with the 6-31G and 6- 
31G** basis sets and the GOOD and TIGHT numerical thresholds 
(see text) controlling numerical precision of the result. The 
lines are fits to the last three and four points, respectively. 



FIG. 5: Total CPU times with increasing order of the response 
for the fifth CPSCF cycle computed as the n + 1 expectation 
value, Eq. II21II . 




number of water molecules 



plemented in the GAMESS quantum chemistry package 
|49| . These static properties have been evaluated at the 
geometries given by Otto et al. 0, and the GAMESS 
results are given to the number of digits provided by that 
program. The MONDOSGF results have been obtained 
both as expectation values, given by Eq. (|21|) . and using 
the non-orthogonal density matrix 2n + 1 rules given in 



Eqs. mm . 

As a benchmark, we have also carried out calculations 
with the linear chain (H20)2o with the VERYTIGHT nu- 
merical thresholding parameters, which employ a 10^^ 
drop tolerance and aim to provide 10 digits of preci- 
sion in the total energy. These VERYTIGHT calculations 
yield = 7.142422 a.u., /J^^^ = -12.033362 a.u. and 
7^^^^ = 1411.425500 a.M.. 
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FIG. 6: Superposition of the magnitudes of the RHF/6-31G 
density matrix derivative elements Dcd, D^dy ^cS ^-nd -Dcd^ 
along the x axis with the separation of basis function centers 
for (H2O)i50. The density matrix derivatives have been con- 
verged to within TIGHT (e.g. a matrix threshold r = 10^^ 
[a.u.]). 
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Linear scaling: 3D water clusters 

Linear scaling computation of the RHF/6-31G and 
RHF/6-31G** second hyperpolarizability, achieved with 
Perturbed Projection, is shown for three-dimensional wa- 
ter clusters in Fig ^ These timings are the total CPU 
time for the fifth CPSCF cycle, including build time 
for J^"''" (ONX and QCTC), iterative construction of 
jyabc (Perturbed Projection via TC2) and all interme- 
diate steps including the congruence transformation. A 
breakdown of the dominant contributions to these totals 
are given in Figs 1 2141 which show timings for Coulomb 
summation (QCTC), Perturbed Projection (TC2), and 
exact exchange (ONX). 

Figure shows the increase in cost associated with 
computing higher order response functions. Correspond- 
ing to this increase in cost, Fig. El shows the magnitude 
of atom-atom blocks of density matrix derivatives up to 
third order as a function of atom-atom distance when 
perturbed by a static electric field. The density response 
shows an approximate exponential decay as a function 
of internuclear distance with the rate of decay slowing 
slightly and the distribution shifted up with increasing 
order in the perturbation. 



DISCUSSION 

In our current formulation, the increase in magnitude 
and reduction of locality in elements of the response func- 



tion makes achieving linear scaling more difficult with 
increasing order in perturbation. Nevertheless, linear 
scaling has been achieved at the HF level of theory up 
to fourth order (i.e. 7) in the total energy for three- 
dimensional systems and non-trivial basis sets. At fourth 
order. Perturbed Projection and exact exchange were 
the dominant costs in solving the CPSCF equations, as 
shown in Figs 01 and 01 For the fourth order Perturbed 
Projection, A^-scaling is achieved between 70 to 110 wa- 
ter molecules, depending on r. Despite a nearly dense 
Y)abc ^ the dominant work in its construction always in- 
volves multiplication with matrices that are significantly 
more sparse, as X°''^'^X^ or X^^Xf. Likewise, A^-scaling 
is achieved between 70 to 90 water molecules for con- 
struction of the Hartree-Fock exchange contribution. In 
this case, the approximate decay of the density matrix 
still leads to linear scaling through ordered skip out lists, 
as described in Ref. [s^. In both cases, the increase in 
response function magnitude is equivalent to tightening 
numerical thresholds, which increases the cost and delays 
the onset of linear scaling. 

In tables IJinl we find that a reduction of the drop tol- 
erance by one order of magnitude leads to an increase in 
precision by 1-2 significant digits, with GOOD and TIGHT 
yielding approximatively 3-4 and 5-6 correct digits inde- 
pendent of the response order. We further observe about 
1 extra digit of accuracy when using the 2?!+ 1 rule. This 
might be expected from the higher order error propaga- 
tion resulting from products of lower order response func- 
tions, relative to evaluation with Eq. (|21|l . which involves 
an error that is always linear in a higher order derivative 
density matrix. 

As shown in Fig. [3 computing the second order re- 
sponse is significantly cheaper than the third order re- 
sponse, and involves an earlier onset of linear scaling. 
Because evaluation of properties with the 2n -f 1 rule is 
of negligible cost relative to solving the CPSCF equa- 
tions, the cost difference for evaluating 7 with the 2n + l 
rule relative to the n + 1 expectation is just the difference 
(roughly 2:3) between the computation of /3 and 7 shown 
in Fig.ini 



CONCLUSIONS 

Linear scaling has been demonstrated for the computa- 
tion of response properties beyond second order in the to- 
tal energy using Perturbed Projection for solution of the 
Coupled-Perturbed Self-Consistent-Field equations. In 
addition, we have provided details of the computational 
method, used three-dimensional systems and non-trivial 
basis sets to demonstrate linear scaling and provided a 
(preliminary) assessment of error control. Perturbed Pro- 
jection for the computation of higher order response func- 
tions is quadratically convergent, simple to implement 
through higher orders and numerically stable. Perturbed 
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Projection is not unique to the Hartree-Fock model, the 
TC2 generator or the MondoSCF iV-scahng algorithms, 
but can be straightforwardly extended to models that in- 
clude exchange correlation (DFT), other other purifica- 
tion schemes such as TRS4 23] as well as other electronic 
structure programs. 

We have shown that response functions (density matrix 
derivatives) through fourth order are local upon global 
electric perturbation, corresponding to an approximate 
exponential decay of matrix elements. However, the mag- 
nitude of the corresponding response functions increases 
with increasing perturbation order, equivalent to tight- 
ening of the matrix drop tolerance, r. While we have not 
attempted to work out a detailed analysis for the prop- 
agation of error, it may be possible to develop a more 
effective thresholding scheme for high orders. In addi- 
tion to being somewhat more accurate, the 2n + 1 rule 
also provides a significantly cheaper alternative to the 
computation of expectation values and an earlier onset 
of linear scaling. 

A similar exponential decay in the first order response 
corresponding to a local nuclear displacement has likewise 
been demonstrated by Ochsenfeld and Head-Gordon p^ . 
This behavior is expected to hold generally for both lo- 
cal and global perturbations to insulating systems. Thus, 
the potential exists for Perturbed Projection to achieve 
linear scaling for a large class of static molecular prop- 
erties within the HF, DFT and hybrid HF/DFT model 
chemistries. Of particular interest, the recently devel- 
oped non-orthogonal density matrixj)erturbation theory 
put forward in a proceeding article 30] may enable linear 
scaling computation of analytic second derivatives, which 
are important in computation of the Hessian matrix. 

This work has been supported by the US Department 
of Energy under contract W-7405-ENG-36 and the ASCI 
project. The Advanced Computing Laboratory of Los 
Alamos National Laboratory is acknowledged. All the 
numerical computations have been performed on com- 
puting resources located at this facility. 
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